library("tidyverse")
library("firatheme")
library("prediction")
library("estimatr")
library("texreg")

# load data
load("~/df_replication.RData")

# Table SI 6
model7 <- lm_robust(click ~ cue*incumbent + issue*incumbent + region + party_name, data = df_replication, se_type="HC2")
texreg(model7, digits = 3, include.ci = FALSE)

# Figure SI 10
a_incumbent <- prediction(model7, at = list(incumbent = 1, cue = "problem_indicators"))
b_incumbent <- prediction(model7, at = list(incumbent = 0, cue = "problem_indicators"))
c_incumbent <- prediction(model7, at = list(incumbent = 1, cue = "public_opinion"))
d_incumbent <- prediction(model7, at = list(incumbent = 0, cue = "public_opinion"))
e_incumbent <- prediction(model7, at = list(incumbent = 1, cue = "rival_parties"))
f_incumbent <- prediction(model7, at = list(incumbent = 0, cue = "rival_parties"))
g_incumbent <- prediction(model7, at = list(incumbent = 1, cue = "the_media"))
h_incumbent <- prediction(model7, at = list(incumbent = 0, cue = "the_media"))

figure10 <- tibble(
  Treatment = c("Problem Indicators","Problem Indicators","Public Opinion", "Public Opinion", "Party Positions", "Party Positions", "Media Reporting", "Media Reporting"),
  Incumbent = c("Yes","No","Yes", "No", "Yes", "No", "Yes", "No"),
  avg = c(mean(a_incumbent$fitted)*100, mean(b_incumbent$fitted)*100, mean(c_incumbent$fitted)*100, mean(d_incumbent$fitted)*100, mean(e_incumbent$fitted)*100, mean(f_incumbent$fitted)*100, mean(g_incumbent$fitted)*100, mean(h_incumbent$fitted)*100), 
  lower = c(mean(a_incumbent$fitted)*100 - mean(a_incumbent$se.fitted)*100, mean(b_incumbent$fitted)*100 - mean(b_incumbent$se.fitted)*100, mean(c_incumbent$fitted)*100 - mean(c_incumbent$se.fitted)*100, mean(d_incumbent$fitted)*100 - mean(d_incumbent$se.fitted)*100, mean(e_incumbent$fitted)*100 - mean(e_incumbent$se.fitted)*100, mean(f_incumbent$fitted)*100 - mean(f_incumbent$se.fitted)*100, mean(g_incumbent$fitted)*100 - mean(g_incumbent$se.fitted)*100, mean(h_incumbent$fitted)*100 - mean(h_incumbent$se.fitted)*100),
  upper = c(mean(a_incumbent$fitted)*100 + mean(a_incumbent$se.fitted)*100, mean(b_incumbent$fitted)*100 + mean(b_incumbent$se.fitted)*100, mean(c_incumbent$fitted)*100 + mean(c_incumbent$se.fitted)*100, mean(d_incumbent$fitted)*100 + mean(d_incumbent$se.fitted)*100, mean(e_incumbent$fitted)*100 + mean(e_incumbent$se.fitted)*100, mean(f_incumbent$fitted)*100 + mean(f_incumbent$se.fitted)*100, mean(g_incumbent$fitted)*100 + mean(g_incumbent$se.fitted)*100, mean(h_incumbent$fitted)*100 + mean(h_incumbent$se.fitted)*100))

ggplot(figure10, aes(x= reorder(Treatment,-avg), y=avg, fill=Incumbent)) +
  geom_bar(stat="identity", position="dodge", color= "black", width=.25) +
  geom_linerange(aes(ymin=lower, ymax=upper), position= position_dodge(.25)) +
  geom_text(aes(y=upper, label=format(round(avg, digits=1), nsmall = 1)), size = 3, position=position_dodge(width=0.7), vjust=-1.75)+
  labs(y="Click Rate (%)", x="", title="", caption = "Note: Error bars represent ± 1 standard error (HC2)", size = 4)+
  scale_fill_manual(values=c("grey80","gray50", "grey80","gray50", "grey30","gray50", "grey80","gray50")) + 
  scale_y_continuous(limits = c(0, 40)) + theme_fira() + ggtitle("a") +
  theme(text = element_text(size=8), legend.position= "bottom", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = 0))


# Figure SI 11
i_incumbent <- prediction(model7, at = list(incumbent = 1, issue = "climate"))
j_incumbent <- prediction(model7, at = list(incumbent = 0, issue = "climate"))
k_incumbent <- prediction(model7, at = list(incumbent = 1, issue = "school"))
l_incumbent <- prediction(model7, at = list(incumbent = 0, issue = "school"))

figure11 <- tibble(
  Problem = c("Climate Change", "Climate Change", "Public Schools","Public Schools"),
  Incumbent = c("Yes","No","Yes", "No"),
  avg = c(mean(i_incumbent$fitted)*100, mean(j_incumbent$fitted)*100, mean(k_incumbent$fitted)*100, mean(l_incumbent$fitted)*100), 
  lower = c(mean(i_incumbent$fitted)*100 - mean(i_incumbent$se.fitted)*100, mean(j_incumbent$fitted)*100 - mean(j_incumbent$se.fitted)*100, mean(k_incumbent$fitted)*100 - mean(k_incumbent$se.fitted)*100, mean(l_incumbent$fitted)*100 - mean(l_incumbent$se.fitted)*100),
  upper = c(mean(i_incumbent$fitted)*100 + mean(i_incumbent$se.fitted)*100, mean(j_incumbent$fitted)*100 + mean(j_incumbent$se.fitted)*100, mean(k_incumbent$fitted)*100 + mean(k_incumbent$se.fitted)*100, mean(l_incumbent$fitted)*100 + mean(l_incumbent$se.fitted)*100))

ggplot(figure11, aes(x= Problem, y=avg, fill=Incumbent)) +
  geom_bar(stat="identity", position="dodge", color= "black", width=.25) +
  geom_linerange(aes(ymin=lower, ymax=upper), position= position_dodge(.25)) +
  geom_text(aes(y=upper, label=format(round(avg, digits=1), nsmall = 1)), size = 3, position=position_dodge(width=0.7), vjust=-1.75)+
  labs(y="Click Rate (%)", x="", title="", caption = "Note: Error bars represent ± 1 standard error (HC2)", size = 4)+
  scale_fill_manual(values=c("grey80","gray50", "grey80","gray50", "grey30","gray50", "grey80","gray50")) + 
  scale_y_continuous(limits = c(0, 40)) + theme_fira() + ggtitle("b") +
  theme(text = element_text(size=8), legend.position= "bottom", plot.caption = element_text(size = 6, color = "gray40", face = "italic"),
        plot.title = element_text(hjust = 0))
